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ABSTRACT 

The statistical study of voids in the matter distribution promises to be an important 
tool for precision cosmology, but there are known discrepancies between theoretical 
models of voids and the voids actually found in large simulations or galaxy surveys. 
The empirical properties of observed voids are also not well understood. In this paper, 
we study voids in an A-body simulation, using the ZOBOV watershed algorithm. As 
in other studies, we use sets of subsampled dark matter particles as tracers to identify 
voids, but we use the full-resolution simulation output to measure dark matter densi¬ 
ties at the identified locations. Voids span a wide range of sizes and densities, but there 
is a clear trend towards larger voids containing deeper density minima, a trend which 
is expected for all watershed void finders. We also find that the tracer density at void 
locations is usually smaller than the true density, and that this relationship depends 
on the sampling density of tracers. We show that fits given in the literature fail to 
match the observed density profiles of voids. The average enclosed density contrast 
within watershed voids varies widely with both the size of the void and the minimum 
density within it, but is always far from the shell crossing threshold expected from 
theoretical models. Voids with deeper density minima also show much broader density 
profiles. We discuss the implications of these results for the excursion set approach to 
modelling such voids. 

Key words: cosmology: observations - large-scale structure of Universe - cosmology: 
theory - methods: data analysis 


1 INTRODUCTION 

The study of large underdense voids in the large-scale matter 
distribution of the Universe has become increasingly impor¬ 
tant in recent years, with the creation of a number of public 
catalogues of voids in galaxy survey data (Pan et al. 2012; 
Sutter et al. 2012; Nadathur & Hotchkiss 2014) and a wide 
variety of statistical analyses based on them. 

Voids are interesting primarily because of the cosmo¬ 
logical information they may contain. Various studies have 
suggested that they could be used to constrain the expan¬ 
sion history of the Universe and the equation of state of 
dark energy (e.g. Ryden 1995; Lee & Park 2009; Bos et al. 
2012; Hamaus et al. 2014), to test modified theories of grav¬ 
ity (Li, Zhao & Koyama 2012; Clampitt, Cai & Li 2013; 
Cai, Padilla & Li 2015; Zivick et al. 2015), to calibrate mea¬ 
surements of galaxy bias (?; Chan, Hamaus & Desjacques 
2014), to constrain initial conditions of structure formation 
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(Kamionkowski, Verde & Jimenez 2009), or to probe more 
exotic theories such as coupled dark energy (Sutter et al. 
2014). The primary void observables used in such studies 
are their abundances and size distributions, the distortion 
of their shapes in redshift space (Alcock & Paczynski 1979), 
their dark matter density profiles, and the void-galaxy or 
void-void position correlations. Given the exciting potential 
applications, a rigorous comparison of theoretical predic¬ 
tions of these properties and those seen for voids in A-body 
simulations and galaxy surveys is very important. 

However, this aim is complicated by the degree of am¬ 
biguity surrounding a very fundamental question: what ex¬ 
actly is a ‘void’? From a theoretical perspective, there is a 
clear answer, provided by the spherical evolution model of 
Sheth & van de Weygaert (2004), who identify voids as those 
non-linear underdense regions which have evolved to reach 
shell crossing. This identification is convenient, as voids can 
then be modelled analogously to collapsed overdense haloes 
using the excursion set formalism (Press & Schechter 1974; 
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Bond et al. 1991; Lacey & Cole 1993), and therefore allows 
clear predictions to be made for void observables. 

In practical terms, however, the definition of a void is 
not so clear. When dealing with either Wbody simulations 
or galaxy survey data, an algorithmic approach is required 
to identify regions as voids, which is complicated by the 
fact that voids are naturally poorly sampled by observable 
tracers, making shot noise a serious issue. A number of dif¬ 
ferent void finders have been proposed (see Colberg et al. 
2008, for a review of methods), which nnfortunately do not 
always agree with each other. Watershed void finders (e.g.. 
Platen, van de Weygaert & Jones 2007; Neyrinck 2008; Sous- 
bie 2011; Cautun, van de Weygaert & Jones 2013) form an 
interesting class of algorithms. They use tessellation tech¬ 
niques (Schaap & van de Weygaert 2000) to reconstruct the 
density field from discrete data points, and the watershed al¬ 
gorithm for creating a void hierarchy. They present a number 
of advantages for practical studies, as they are more resilient 
to shot noise in the density reconstruction (though see also 
e.g. Elyiv et al. 2015; Shandarin & Medvedev 2014, for other 
interesting proposals), and do not make prior assumptions 
about void geometries. They are also the most commonly 
used. Watershed voids may therefore be considered a rea¬ 
sonable practical definition of a void. 

However, watershed algorithms make no reference to 
shell crossing, which is the defining characteric of theoretical 
models. The obvious question is therefore how, or whether, 
these two void definitions are related to each other. The an¬ 
swer is important to the practical use of watershed voids in 
cosmology, as well as to the development of further theo¬ 
retical predictions. A number of studies (e.g., Sutter et al. 
2014; Chan et al. 2014; Pisani et al. 2014; Chongchltnan 
2015; Achitouv, Neyrinck & Paranjape 2015) which apply 
the Sheth & van de Weygaert (2004) formalism to describe 
watershed voids assume that the two approaches describe 
the same or closely related objects. There is known to be a 
significant disagreement between model predictions for void 
abundances as a function of their size and results obtained 
for watershed voids in dark matter simulations. This can 
be partially resolved (at least at large void sizes) by the 
ad hoc assumption that shell crossing and void formation 
occur at less extreme densities than predicted by the spher¬ 
ical model, but as we will also show, such an approach lacks 
self-consistency. A more direct comparison of void properties 
with the model is therefore desirable. 

At the same time, a number of properties of watershed 
voids remain imperfectly understood. A generic property of 
watershed void finders is that voids containing the deepest 
density minima should have the largest sizes. Yet the fits 
provided by Hamaus, Sutter & Wandelt (2014) to describe 
density profiles about void centres appears to suggest the op¬ 
posite behaviour (though note that the applicability of this 
fitting form is not universally accepted, e.g. Nadathur et al. 
2015). Perhaps related to this problem is the question of how 
to define the ‘centre’ of a void — which is also important for 
void correlation studies. The standard procedure assigns the 
centre to a weighted average of the positions of the tracers 
of mass within a void (e.g. Lavaux & Wandelt 2012; Sutter 
et al. 2012; Nadathur & Hotchkiss 2014; Sutter et al. 2014), 
but it is not clear that this will correctly identify the re¬ 
gion with the greatest absence of mass. Another interesting 
question is how densities reconstructed from discrete tracer 


distributions relate to the true underlying density field. In 
studies of voids from simulations, the full simulation output 
is typically randomly down-sampled to provide a set of trac¬ 
ers, which can have a dramatic effect on the recovered void 
properties (Sutter et al. 2014). 

Our goals in this paper are two-fold. We wish to under¬ 
stand the relationship between watershed voids and theoret¬ 
ical models. To do this, we move beyond the fitting of void 
number functions alone and identify other important char¬ 
acteristics of voids which can be used to test the assumption 
of shell crossing more broadly. We also want to empirically 
examine the properties of watershed voids in simulations in 
order to understand the working of the algorithm and clarify 
some of the issues above. 

To do so we make use of the popular ZOBOV watershed 
algorithm (Neyrinck 2008). To enable comparison of our re¬ 
sults with others in the literature, we will mostly use the 
options for ZOBOV implemented in the VIDE toolkit (Sut¬ 
ter et al. 2015). We analyse voids identified using randomly 
subsampled dark matter particles as tracers, and relate them 
to true densities determined from the full resolution simula¬ 
tion output. We propose a new definition of the void centre, 
which is designed to better identify the true location of the 
underdensity within the void. The details of the Y-body 
simulation, the watershed algorithm and the methods for 
identifying void centres and measuring density profiles are 
described in Section 2. 

Section 3 provides a summary of the spherical model 
for void evolution, which we use to extract general identify¬ 
ing characteristics of shell-crossed voids for comparison with 
simulation results. This comparison is performed in Sec¬ 
tion 4, where we also outline the general properties of voids 
in our simulation. We show that larger voids do correspond 
to deeper density minima, as expected for the watershed 
algorithm. We examine the viability of the fitting formula 
of Hamaus et al. (2014) to describe the density profiles of 
voids, and show that the fits described in that paper do not 
provide a good quantitative or qualitative description of the 
variation of the average profile within the void population. 
In addition, subsampled tracers almost always overestimate 
the density contrast in voids. All of these results have prac¬ 
tical implications for future studies that use watershed void 
finders. We compare these results from simulated voids to 
theory and argue that there is no evidence that watershed 
void finders in general, and VIDE and ZOBOV in particular, 
satisfy the primary defining criteria of the Sheth & van de 
Weygaert (2004) model. This leads us to reassess the viabil¬ 
ity of describing watershed voids using existing theoretical 
techniques. We summarize and conclude in Section 5. 

2 NUMERICAL METHODS 
2.1 Simulation 

In this paper we use Y-body simulation data from 
the MultiDark Runl (MDRl) release of the MultiDark 
project (Prada et al. 2012).^ MDRl uses an Adaptive- 
Refinement-Tree (ART) code, based on adaptive mesh re¬ 
finement, to simulate 2048^ dark matter particles within 

^ Publicly available from www.cosmoslm.org. 
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a cubic volume of 1 {h~^Gpc)^, in a ACDM cosmo¬ 
logical model with parameters (Jim, JIa, fib, h,, ris, ug) = 
(0.27,0.73,0.0469,0.7,0.95,0.82). The simulation has mass 
resolution nip = 8.721 x 10® h~^M q and force resolution 
7 /i“^kpc. Initial conditions were set using the Zeldovich 
approximation at redshift z = 65. 

From the full particle output at redshift 0 we randomly 
subsample the dark matter particles down to a nnmber den¬ 
sity of n = 3.2 X 10“® h®Mpc“®, similar to that of typical 
galaxy samples (e.g. Zehavi et al. 2011). This corresponds to 
a mean nearest-neighbonr separation of ~ 7 /i“^Mpc. 

We refer to the resnlting sample as the Main sample, and use 
these particles as tracers for the void finding. In addition, 
we have used a control sample with a higher tracer density 
2 X 10“® /i®Mpc“®, which we refer to as the Dense sample. 
However, our primary conclusions regarding the properties 
of watershed voids do not depend strongly on the tracer 
number density. Therefore unless otherwise stated, all re¬ 
sults presented in this paper refer to voids from the Main 
sample. 

Note that a random subsampling of tracers introduces 
shot noise but does not change the fundamental clustering 
properties of the dark matter field. Therefore despite having 
the same tracer number density, the properties of voids in 
such subsampled tracers and those in the galaxy distribution 
would not be expected to be (and are not) the same, since 
galaxies are biased tracers of the matter density. However 
for our purposes of understanding the general properties of 
watershed voids in this work, a random subsampling is suf¬ 
ficient. We consider the effects of galaxy bias separately in 
a companion paper (Nadathur & Hotchkiss 2015). 

Although our tracers are themselves dark matter parti¬ 
cles, as the subsampling procedure increases shot noise we 
will distinguish between the tracer number density, and the 
underlying dark matter density. The dark matter density in 
MDRl is determined from the full resolution particle out¬ 
put of the simulation at redshift 0, by using a cloud-in-cell 
interpolation on a 1024® grid, followed by smoothing with a 
Gaussian kernel with width equal to one grid cell. The sub- 
Mpc resolution of this grid is much smaller than the typical 
void size scales, so that this procedure in effect provides 
a continuous underlying dark matter density, of which the 
subsampled tracer particles are an approximately Poisson 
realization. 

In the following, we will reserve the symbols p and 
A for dark matter densities determined using this gridded 
smoothed density field, and use the symbols n and A„ for 
the equivalent quantities determined from the tracer number 
densities. 


2.2 Void finding 

To identify voids in the dark matter particle distribution, 
we make use of the ZOBOV watershed void finder (Neyrinck 
2008), with the options implemented in the VIDE toolkit 
(Sutter et al. 2015). Although there are known issues with 
the application of VIDE to galaxy survey data with irregular 
survey volumes and masks (Nadathur & Hotchkiss 2014), 
when dealing with a simulation cube with periodic boundary 
conditions these do not present a problem. However, note 
that in some cases, especially the definition of the void centre 
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described below, we use our own modification of the ZOBOV 
algorithm. 

ZOBOV works by reconstructing the density field based 
on a Voronoi tessellation of the simulation cube around the 
discrete distribution of tracer particles. The tessellation as¬ 
sociates each particle with a Voronoi cell consisting of the 
region of space closer to it than to any other particle. The 
volume of the Voronoi cell i relative to the mean volume 
is then used to estimate the local tracer number density 
Ui at the particle location. This reconstruction, known as 
the Voronoi tessellation field estimator (VTFE), is naturally 
scale-adaptive and thus far more resilient against shot noise 
effects than naive counts-in-cells measurements. 

After reconstructing the density field, the algorithm 
identifies all local minima of the reconstructed density field 
and determines the “catchment basins” around each min¬ 
imum, known as zones. Zones are then merged to form a 
nested hierarchy of voids according to the watershed princi¬ 
ple (Platen et al. 2007), such that the zone with the smallest 
minimum density rimin then acquires neighbouring higher- 
density zones as sub-voids, in increasing order of the mini¬ 
mum density on the watershed ridge separating them. For 
each void thus formed, we define an effective void radius 
as the radius of a sphere with equivalent volume V, 

/3 

• ( 1 ) 

Even in the absence of any merging, deeper density 
minima typically correspond to zones of larger volume and 
thus larger Rv However, the watershed merging procedure 
also ensures that voids with deepest density minima contain 
greater numbers of merged sub-voids and therefore have the 
largest sizes. This correlation of minimum density and void 
size is a common property of all watershed void-finders and 
is not unique to the ZOBOV algorithm. 

To avoid excessive merging leading to essentially infinite 
void sizes, VIDE imposes a restriction preventing the merger 
of two zones unless the minimum link density along the wa¬ 
tershed ridge separating them satisfies nunk < n-max, where 
Umax = 0.2n. This condition applies only to the lowest den¬ 
sity point on such a ridge, and does not prevent voids from 
containing regions of much higher densities. The value of 0.2 
therefore has no theoretical motivation, and rimax should be 
considered an arbitrary free parameter. Indeed alternative 
values of Umax have been considered in other works (Na¬ 
dathur & Hotchkiss 2014; Hotchkiss et al. 2015; Achitouv 
et al. 2015; Nadathur et al. 2015), and properties such as 
the abundance of root-level voids, the distribution of void 
sizes and void density profiles will all depend on the value 
chosen. A fuller discussion of the effects of this arbitrary 
choice is provided by Nadathur & Hotchkiss (2015). How¬ 
ever, for ease of comparison with previous results we shall 
restrict ourselves in this paper to the default value hard¬ 
coded in VIDE. 

Further selection cuts might be desirable at this stage, 
since ZOBOV simply reports all local density minima as po¬ 
tential voids, without regard to the value of the minimum 
density within them or any reference to shell crossing. VIDE 
provides an optional selection cut which purports to remove 
those voids which have a tracer number density within a de¬ 
fined central region greater than 0.2 times the mean. How¬ 
ever, this density is measured by naive number counting on a 
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scale smaller than the mean inter-particle separation. There¬ 
fore, as pointed out by Nadathur & Hotchkiss (2013), it is 
very badly affected by shot noise and the values determined 
by VIDE are almost completely uncorrelated with the true 
central density. In any case this cut only excludes a small 
fraction of final voids, so we do not apply it. We also do not 
apply the much more conservative cuts on rimin suggested by 
Nadathur & Hotchkiss (2014), Hotchkiss et al. (2015) and 
Nadathur et al. (2015). 

A selection cut based on the void radius has some¬ 
times been advocated in the literature, to remove voids with 
Rv < which are claimed to be below the resolution 

limit. In fact the adaptive nature of the tessellation means 
that ZOBOV automatically excludes small voids below its 
resolution limit, as we show in Section 4. Therefore no fur¬ 
ther cut on Rv is necessary. 

By applying these criteria we find in total 27 450 voids 
in the Main tracer sample. Of these, 26 919 are root-level 
voids in the hierarchy, i.e. they are not subvoids of any par¬ 
ent voids and their volumes do not overlap each other. 


tracer density rimin), and examine the tessellation output to 
identify all Voronoi cells adjacent to it. From this set we 
select the lowest density neighbouring particle, and then, in 
order of increasing density, two other particles that are ad¬ 
jacent to both the core particle and the previous selections. 
This provides us with the four lowest density mutually adja¬ 
cent Voronoi cells in the void; we now dehne the void centre 
to lie at the point of intersection of these four cells, which is 
also the circumcentre of the tetrahedron formed by the four 
tracer particles. This point represents the location within 
the void that is maximally distant from all tracers. 

We shall refer to this alternative definition of the void 
centre as the circumcentre and denote its location by In 
Section 4 we show that both the number density of tracers 
and the underlying dark matter density are indeed signif¬ 
icantly lower at the circumcentre than the barycentre. In 
Appendix B we also show that the circumcentre location is 
more resilient to shot noise effects arising due to subsam- 
pling. 


2.3 Void centres 


Since voids obtained from the watershed algorithm have ar¬ 
bitrary shapes, different prescriptions may be used to define 
the location of the void ‘centre’. The most commonly used 
definition, which is also the definition implemented in VIDE, 
is the volume-weighted barycentre of the void member par¬ 
ticles 


X 


be 

V 



( 2 ) 


where Xi is the position of the ith particle, Vi is the volume 
of it’s corresponding Voronoi cell, and the sum runs over all 
member particles of void v. 

In the low-density interior of a void, Voronoi cells in the 
tessellation are typically greatly elongated, and the particles 
contained within them lie far from their geometrical centres. 
This means that the position Xi corresponding to each cell is 
an imprecise measure of the location of the cell. In addition, 
watershed voids contain a great number of member parti¬ 
cles - the median number for our void sample is 82, and 
many voids contain several hundreds - the vast majority 
of which reside in the overdense walls and filaments on the 
outskirts of voids. A combination of these two factors means 
that although the barycentre position defined by equation 2 
is roughly symmetrically located with respect to the over- 
dense void walls, it is typically very far from the position of 
minimum density. This is because the barycentre dehnition 
is fundamentally based on the locations where tracers are 
present, rather than locations from where they are absent. 
A consequence of this is that the location of the barycentre 
is very sensitive to the sub-void fraction and thus to the ar¬ 
bitrary condition controlling void merging described above. 

For some purposes, it may be more logical to dehne the 
void centre to coincide with the location of the minimum 
density within it. This is particularly important when void 
centre locations are subsequently used to measure density- 
dependent effects, such as void leasing or ISW contribu¬ 
tions. To achieve this, we adopt the following procedure. 
We identify the core particle of the void as the particle with 
the largest Voronoi cell (i.e., corresponding to the minimum 


2.4 Density profile determination 


A fundamental quantity of interest is the average distribu¬ 
tion of tracers and dark matter about the void centre, and 
the variation of this distribution with void properties. We 
study this behaviour by constructing stacked density pro- 
Hles for subsets of voids satisfying different criteria. To do 
so we rescale distances within each qualifying void in units 
of the void radius Rv, and then estimate the average den¬ 
sity in the stack in concentric spherical shells about the void 
centre. 

Estimating the tracer number density in this way is 
complicated by shot noise effects, since the interiors of voids 
by dehnition contain very few tracer particles which can 
be used for number density measurements. Nadathur et al. 
(2015) showed that an unbiased estimate accounting for 
Poisson noise can be obtained using the volume-weighted 
estimator for the average number density in the jth radial 
shell. 




( 3 ) 


where the jth shell has width Af in units of the rescaled 
radial distance f for each void, V^’ is the true volume of the 
jth shell of the ith void and V/ is the number of tracer par¬ 
ticles contained within it, and the sum over i runs over all 
voids included in the stack. Note that the in this formula 
includes all tracer particles within the shell, not just those 
that are identihed as members of the void by the water¬ 
shed algorithm. Under the assumption that the individual 
numbers are Poisson realizations of the true underly¬ 
ing density, the error in equation 3 can then be estimated 
at any desired confidence level directly from the definition 
of the Poisson distribution. In this paper plotted errorbars 
indicate the 68% conhdence limits on n. 

As the resolution of the dark matter density held is 
much hner than the typical void size, no such complications 
are required when estimating p over the stack of voids. We 
simply sample the dark matter density at all grid points 
contained within the radial shell and calculate the mean 
and standard deviation of the values obtained. We present 
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our results for the stacked dark matter densities in terms 
of the average total enclosed density within a radius r, 1 + 
A(r), since this allows a more direct contact with the theory 
described in Section 3. 

All density proHles are calculated out to three times the 
void radius, and are measured in radial bin steps of 0.1 times 
the radius. 


3 EXCURSION SET MODELS OE VOIDS 


Most existing theoretical descriptions of voids are derived 
from the framework presented by Sheth & van de Weygaert 
(2004). This in turn derives from the original excursion set 
approach of Press & Schechter (1974); Epstein (1983); Bond 
et al. (1991) and is based on the model of spherical evolution 
of mass shells (Gunn & Gott 1972; Lilje & Lahav 1991). In 
this Section we briefly summarize such models in order to 
highlight the key areas of comparison with the results of 
watershed void finders. 

In this picture the evolution of a spherical mass shell of 
radius r is determined by the total enclosed density contrast 
within the radius of the shell at time t, 1 + A(r, t), where 


A{r,t) 


rr 

' p{y,t) f 

Jo 

[ m \ 


y^dy, 


( 4 ) 


and by the time evolution of the cosmological density pa¬ 
rameter fl(t). Underdense spherical regions contain a density 
deficit (i.e., A{r,t) < 0) which causes shells to expand out¬ 
wards. This deficit is stronger for inner shells, which there¬ 
fore expand faster than outer shells, and mass evacuated 
from the centre of the underdensity begins to pile up at its 
edges. For a steep enough starting density profile, at some 
point in the evolution inner shells catch up with shells which 
were initially further out from them, in an event known as 
shell crossing. The moment of shell crossing marks a transi¬ 
tion in the evolution of the underdensity, as it subsequently 
expands outwards self-similarly (Suto, Sato & Sato 1984; 
Fillmore & Goldreich 1984; Bertschinger 1985). 

Within the spherical model, it can be shown that shell 
crossing occurs when the average density enclosed within 
the void is 


Penc/p = 1 -I- A(r, t) ~ 0.2. (5) 

This corresponds to a linearly extrapolated average density 
contrast of 

Alin = Jv -2.71 , (6) 

at the epoch of shell crossing, with this value independent of 
radius r and largely independent of the cosmological param¬ 
eters governing the background evolution. This is analogous 
to the case of spherical collapse of clusters, which occurs 
above a linear overdensity threshold of Aun = 5c — 1.69. 

Following Blumenthal et al. (1992); Dubinski et al. 
(1993), Sheth & van de Weygaert (2004) then identify the 
population of voids with only those mature evolved under¬ 
densities that have reached the stage of shell crossing. If the 
initial Gaussian density fluctuation held is smoothed on a 
range of different smoothing scales R, this physical picture 
identihes huctuations which exceed the density threshold 5v 
on smoothing scale R with potential voids of radius R to¬ 
day. These huctuations can be characterized by their depth 
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in units of the rms huctuation of the density held on scale 
Rf 

v = &l/al{R), (7) 


where (Jo{R) is one of the set of spectral moments 
r 

a^{R) = j -^W\kR)P[k)dk , ( 8 ) 

with P{k) the power spectrum of the unsmoothed density 
huctuation held and W{kR) the smoothing hlter. 

At this point, the abundance and size distribution of 
voids can be predicted by a number of models of varying 
degrees of sophistication. Sheth & van de Weygaert (2004) 
use the excursion set approach (e.g.. Bond et al. 1991; Sheth 
1998) to account for huctuations which cross the 5v thresh¬ 
old on some small scale but are overdense with Ann > 5c on 
some larger scale. Such underdensities would be crushed by 
the collapse of the surrounding cluster and so would not 
be visible as voids today (the void-in-cloud effect). This 
amounts to a two-barrier problem. According to this model, 
assuming void number density is conserved on evolving from 
Lagrangian to Eulerian space, this number density can be 
expressed as a function of the Eulerian void radius Rv (e.g., 
Jennings, Li & Hu 2013; Chan et al. 2014) as 



Here the Lagrangian radius Rl = 0.58Rv, a relationship 
determined by the shell crossing condition above. 

However, it is well known that equation 9 does not pro¬ 
vide a good ht to the distribution of voids found by wa¬ 
tershed algorithms, since it predicts a sharp cutoff in void 
sizes above ~ 5 /i“^Mpc, much smaller than observed for 
watershed voids. A number of studies (Jennings et al. 2013; 
Sutter et al. 2014; Chan et al. 2014; Pisani et al. 2015)) 
have attempted to improve fits by relaxing the shell cross¬ 
ing condition 5v = —2.71 and treating 5v as a free parame¬ 
ter instead.^ This procedure is not justified by any specific 
theoretical model. However, if 5v is allowed to vary, self- 
consistency requires that the relationship between Eulerian 
and Lagrangian radius be correspondingly modified to 


Rl 



(1 - 5v/c)=/3 ’ 


( 12 ) 


where c 1.594 (Bernardeau 1994; Jennings et al. 2013). 
Despite this modification, equation 9 with a variable still 
fails to describe the distribution of small voids, and the fit 


^ Note that the definition u = Sv/ao is also commonly used. In 
this case equation 10 would need to be appropriately modified, 
as done by Chan et al. (2014). 

^ Jennings et al. (2013) also propose an alternative adaptation of 
this model, but this cuts off the distribution at even smaller Ry, 
so would make the discrepancy worse. 
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R, [/i-^Mpc] 


Figure 1. The differential number density of voids in simulation 
as a function of their size, for both simulation samples. Error 
bars are calculated assuming the void numbers in each bin are 
Poisson distributed. The dashed line shows the best fit of the 
Sheth & van de Weygaert (2004) model to the > 25 /r“^Mpc 
data, with Sv = —0.40. The solid line shows an exponential cutoff 
model which describes the same data better. 


values of Sv for large voids vary widely. Chan et al. (2014) 
obtain — — 1 and find little redshift dependence of this 
value, contrary to theoretical expectation — however, they 
keep Rl = 0.58i?„ fixed when varying 5^ instead of using 
equation 12, so their model is not self-consistent. On the 
other hand, Sutter et al. (2014) find a range of different Sv 
values for voids from different samples, ranging from —0.26 
to —0.5. Pisani et al. (2015) quote Sv = —0.45. It is hard to 
conceive of an explanation for why Sv should lie so far from 
the theoretical prediction if the shell crossing model were 
true. 

In any case, insofar as allowing Sv to vary allows a fit 
to the distribution of the largest voids, it does so simply 
by replicating an exponential cutoff in the distribution at 
large Rv To demonstrate this, in Fig. 1 we show the dis¬ 
tribution of void sizes obtained from both our simulation 
samples, together with the prediction obtained from eqs. 9 
and 12 with value Sv = —0.40, and a simple exponential 
curve oc exp(—7?S ®°). These numerical values represent the 
best fits obtained from fitting to the Rv > 25 /i“^Mpc data 
from the higher resolution Dense sample (note that the size 
distribution at large Rv is itself resolution-dependent). The 
minimum radius cut is imposed since neither model can fit 
the data at all scales. The best-fit parameters are sensitive to 
the exact choice of this cut, but the relative quality of the 
fits does not change significantly. Despite having an extra 
parameter, the exponential cutoff model significantly out¬ 
performs the modified excursion set model on the basis of 
the Akaike Information Criterion (Akaike 1974). 

An approximate exponential cutoff at large void radii is 
a rather generic feature of alternative descriptions of voids. 


and would also apply if for instance voids were modelled 
simply as minima of a Gaussian density field at a fixed 
smoothing scale (Bardeen et al. 1986; Flender, Hotchkiss & 
Nadathur 2013; Nadathur et al. 2014). We do not intend to 
attempt a fully-fledged alternative theoretical description of 
voids here; instead our point is that the ability or otherwise 
of equation 9 with variable (5v to fit the void distribution in 
a limited size range should not be taken as evidence that 
the excursion set model provides a good description of wa¬ 
tershed voids. 

However, other aspects of the excursion set model can 
also be directly tested. The crucial ingredient of the model 
is identification of underdensities as voids only if they have 
undergone shell crossing. Various modifications of the model 
(e.g., Paranjape, Lam & Sheth 2012; Musso & Sheth 2012; 
Paranjape & Sheth 2012; Jennings et al. 2013; Achitouv 
et al. 2015) do not change this fundamental picture. Since 
watershed algorithms in general make no explicit reference 
to shell crossing when defining a ‘void’, it is desirable to test 
this assumption much more directly than through the ad 
hoc fitting of the number function described above. 

A key property of shell crossing is that it occurs at the 
same enclosed density contrast A for all voids, irrespective of 
their size. (This is also what justifies the use of a single Sv in 
fitting void abundances.) Adding additional complexities to 
the spherical evolution model, such as the effect of a shear 
field, could relax the condition A = —0.8 to some extent, 
and may introduce a small scatter in the values of A over 
the void population. Nevertheless, strong variation of the 
enclosed density with properties of watershed voids would 
be a clear sign that they do not correspond to similar shell- 
crossed objects. 

A related property of shell-crossed voids is that if the 
enclosed density contrast is to be the same for voids of all 
sizes, equation 7 requires that larger voids must correspond 
to more extreme fluctuations of the parameter v. It can 
be shown that this in turn means that larger voids should 
on average correspond to shallower but broader initial den¬ 
sity profiles Sir), while smaller voids correspond to deeper 
and steeper profiles. That is, smaller shell-crossed voids 
should contain deeper density minima than large voids. ^ The 
analogous situation for collapsing haloes is that the most 
massive haloes should be the least centrally concentrated, 
which is indeed the case (e.g. Navarro, Frenk & White 1996, 
1997). More generally, voids with the deepest density min¬ 
ima should have the steepest density profiles, and vice versa. 

These qualitative properties provide clear tests of the 
assumption that watershed voids have undergone shell cross¬ 
ing. However, as we show in the next Section, neither of them 
hold true for the voids obtained using VIDE and ZOBOV, nor 
should one expect them to hold for other watershed void 
finders. 

4 PROPERTIES OF WATERSHED VOIDS 
4.1 Sizes and densities 

Fig. 2 shows the distribution of void sizes and minimum 
tracer number densities for all voids in our Main tracer sam- 

^ We thank Ravi Sheth for drawing our attention to this point. 
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Figure 2. The distribution of the minimum tracer number densities within voids and void sizes in the Main sample. There is a clear 
trend towards increasing void size as the minimum density decreases. The dotted lines show the contours enclosing 95% and 99% of all 
‘voids’ identified in a random uniform distribution of points with the same number density and in the same volume. The arrow indicates 
the value Rv = roughly the mean inter-particle separation, which has sometimes been suggested as minimum size cut. The dashed 

line shows the true minimum achievable void size resolution as a function nmin: most voids automatically lie well away from this limit. 


pie. It is immediately obvious that lower minimum number 
densities are correlated with larger void sizes, as we argued 
would always be the case for watershed void finders in gen¬ 
eral and ZOBOV in particular. The characteristic banana¬ 
shaped distribution is similar to that found by Nadathur 
et al. (2015); Nadathur & Hotchkiss (2015), indicating that 
this is a universal property of the void finder, independent of 
whether the tracer particles used are dark matter particles 
or galaxies in haloes. It is also noteworthy that neither VIDE 
nor ZOBOV impose any restriction on the allowed minimum 
densities, resulting in a range of rimin values extending all 
the way up to the mean. 

The minimum achievable void size resolution is dictated 
by the process of reconstructing the density field from the 
Voronoi tessellation described in Section 2.2, which sets a 
natural cutoff of Rv,min = (3/47r)^^^ vn (rimin/n)”^^^, where 
rjv = ~ 7 /i“^Mpc is roughly the mean inter-particle 

separation. This cutoff is shown by the dashed black line, 
and rjv by the vertical arrow. In fact most voids naturally 
lie well away from this limit, simply because most zones 
contain several tracer particles. This means that the selec¬ 
tion criterion advocated by some studies (e.g. Sut¬ 

ter et al. 2012; Sutter et al. 2014) has no practical effect, 
whereas a tighter criterion R^ > 2rN (Hamaus et al. 2014) 
is unnecessarily conservative. 

Also shown are contours showing the 95 and 99 per 
cent confidence limit contours for the distribution of spu¬ 
rious ‘voids’ identified by the same algorithm in a random 
uniform distribution of points with the same volume and 
same number density as the Main sample. There is clearly a 
considerable overlap between the two distributions, but care 
is required in its interpretation, as P ((nmin, i7„)|Poisson) is 
not the same as P (Poisson|(nmin, Pv))- In the absence of 
information on the true dark matter content of such voids. 


conservative cuts to the void catalogue based on the prop¬ 
erties of Poisson voids have previously been advocated (e.g. 
Neyrinck 2008; Nadathur & Hotchkiss 2014; Hotchkiss et al. 
2015), but these may use available data sub-optimally. Since 
we have access to this information from the simulation, 
we analyse all voids without imposing such cuts a priori, 
and in fact we find that while Poisson contamination in¬ 
creases in the overlap region, statistically speaking voids of 
all (rimin, Rv) values on average correspond to true dark mat¬ 
ter underdensities. Similarly, we find that low values of the 
density ratio r (Neyrinck 2008) also do not serve as a reliable 
indicator of Poisson contamination. This issue is discussed 
further in Appendix B. 

Another noteworthy aspect of Fig. 2 is the apparent 
saturation of the minimum densities within voids with in¬ 
creasing Rv This is a consequence of the finite tracer num¬ 
ber density, and the saturation value is dependent on the 
mean density n. Subsampling tracer particles lowers n and 
thus reduces the apparent tracer density contrast in voids. 
Conversely, Nadathur & Hotchkiss (2015) show that at the 
same mean tracer density, more highly biased tracers result 
in much lower values of rimin within voids. 

Given the uncertainties associated with the tracer num¬ 
ber density discussed below, the true dark matter density 
at void locations is perhaps a more informative quantity. To 
measure this we make use of the dark matter density field de¬ 
scribed in 2.1 and simply measure its value p in the grid cell 
corresponding to the position of the void centre. Fig. 3 shows 
the binned average central densities as a function of the void 
radius Rv , for both Dense and Main voids, and for both def¬ 
initions of the void centre described in Section 2.3. The error 
bars represent the 2a uncertainty in the mean. As expected, 
in all cases the circumcentre definition is a superior indicator 
of the location of minimum density within the void. It is also 
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Figure 4. The distribution of the average tracer number densities within voids and void sizes in the Main sample. The dotted lines 
show the contours enclosing 95% and 99% of voids in random distributions as in Fig. 2, and the arrow indicates the approximate mean 
inter-particle separation. 
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Figure 3. Binned average values of the dark matter density at 
the location of the void centre, as a function of the void size. The 
top and bottom panels show data for voids in the Dense and Main 
samples respectively. Circles (green) and squares (red) refer to the 
two alternative definitions of the void centre; the circumcentre is 
clearly a better locator of the true minimum density in the void. 


clear that this minimum density decreases with increasing 
void size, contrary to the excursion set prediction. Curiously, 
for voids in the Dense sample, the density at the barycentre 
first decreases and then increases with This is because 
in this case voids with large Ru tend to be formed from the 
merger of several sub-voids. As such sub-voids necessarily 
correspond to shallower density minima they do not affect 
the location of the minimum density, but they do shift the 


location of the barycentre via equation 2. In contrast, the 
location of the circumcentre is independent of the sub-void 
fraction and the choice of criteria to control void merging. 

Fig. 4 shows the distribution of the average tracer den¬ 
sity Uavg within the void, calculated from the number of 
void member particles and the void volume, and R^. As 
pointed out by Nadathur & Hotchkiss (2014); Achitouv et al. 
(2015), riavg is typically > 1 and much larger than Umin, 
simply because the watershed definition means that voids 
always extend to include high density regions on the sepa¬ 
rating ridges. This feature of ZOBOV and VIDE indirectly 
demonstrates that most tracers in identihed voids reside in 
overdensities, which explains why the barycentre is a poor 
locator of the minimum underdensity, and also suggests that 
when riavg > 1 the void radius is a signihcant overestimate 
of the size of the true underdense region. 

It is worth noting that a selection cut on the minimum 
void radius alone is a sub-optimal way of excluding voids 
that are on average overdense, since it would eliminate many 
with the lowest riavg values as well. 


4.2 Tracer density versus dark matter density 

The relationship between the tracer number density and 
the true underlying dark matter density in the simulation 
is also of interest. Even though the tracers in our case are 
down-sampled dark matter particles, we find that these two 
quantities are in general not the same. There is already an 
inherent shot noise in the number densities of dark matter 
particles arising from the fact that they constitute a discrete 
realization of the underlying continuous density field. Ran¬ 
domly down-sampling the dark matter particles enhances 
this shot noise, meaning that, particularly in void regions, 
tracer number densities tend to be larger than the true dark 
matter density. This problem is to a large extent mitigated 
by the self-adaptive nature of the Voronoi tessellation, but 
cannot be completely removed. A second consequence of re- 
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Figure 5. Binned average values of the true dark matter density 
at the location of the void centre, as a function of the minimum 
tracer number density within the void, for both definitions of the 
void centre. The top and bottom panels show the data for voids 
from the Dense and Main tracer samples respectively. A 45° line 
is shown for reference in each case. 


ducing the total number of tracers is to increase the effec¬ 
tive smoothing scale at which tracer number densities are 
measured. As mentioned above, for the VTFE density re¬ 
construction this scale is ~ rjv ~ 7 h“^Mpc, whereas the 
dark matter density field is smoothed at a scale of < lh~^ 
Mpc. This change of smoothing scales acts in the opposite 
direction, tending to make dark matter underdensities ap¬ 
pear shallower in riniin. 

The relative strengths of these two effects are illustrated 
in Fig. 5, which shows the relationship between the dark 
matter density p at the position of the void centre and the 
minimum tracer number density within the void as deter¬ 
mined from the tessellation, for both dehnitions of the void 
centre, and for voids from both the Main and Dense tracer 
samples. That p exceeds rimin at the barycentre for both 
samples is to be expected, since the barycentre typically lies 
quite far from the location of the tracer density minimum. 
But even at the circumcentre, which is guaranteed to lie in 
the region of minimum tracer number density, the tracer 
number density does not accurately reflect the dark matter 
density, although there is a clear linear relationship between 
p and rimin- Particularly in the most underdense voids, shot 
noise enhancement dominates, causing rimin/fi < pmin/p- 
However, for the shallowest voids the smoothing effect be¬ 
comes more important, reversing the relationship. The same 
qualitative trends are seen for both the Dense and Main 
tracer samples, and more generally hold for any tracer pop¬ 
ulation obtained from subsampling dark matter particles. 


4.3 Density profiles 

We now turn to the distribution of tracer particles and dark 
matter around void centres. Anticipating that the form of 
the density profiles will depend on the void size, we hrst ex¬ 
amine the average profiles for stacks of voids within different 
ranges of i?„, chosen such that each stack contains an equal 
number of voids. 

Fig. 6 shows the resulting tracer number density profiles 
for stacks centred on the void barycentres and circumcen- 
tres in the left and right panel respectively. The barycentre 
stacks show a strong trend for decreasing central density 
as the void size increases, as is expected from Fig. 2. Voids 
are generally surrounded by overdense walls, which are much 
higher for small voids than for large ones. The general asym¬ 
metry of the circumcentre location with respect to particles 
in the void walls is also apparent in the fact that the stacked 
profiles about this location are less able to resolve the high 
densities in these walls. On the other hand, central den¬ 
sities are much lower for the circumcentre stacks. This is 
the essential tradeoff between the two centre definitions: the 
barycentre has a greater degree of symmetry with respect 
to the surrounding overdensities, whereas the circumcentre 
identihes the true location of the iinderdensity. 

A curious feature is apparent in the stacked barycen¬ 
tre profiles at smaller void radii: the tracer density does 
not show a minimum at the void barycentre, but instead at 
some distance away from the centre. In fact for the smallest 
voids the average tracer density at the barycentre is indis¬ 
tinguishable from the mean. As argued in Section 2.3, the 
void barycentre is always displaced away from the location 
of the minimum tracer density within the void; in particu¬ 
lar, for small voids the barycentre is often at or very close to 
the location of a tracer particle within the void. The tracer 
number density n(r) is measured by naively counting the 
numbers of tracer particles within volumes on scales gener¬ 
ally much smaller than the mean inter-particle separation. 
As a result, when the barycentre location is close to a tracer 
particle, high central values for n(r) are obtained. 

The converse effect can be seen by comparison with the 
right panel of Fig. 6, which shows the stacked density pro¬ 
files for the same voids, but based around the void circum¬ 
centre. The circumcentre is also a special point, as it is by 
construction as far as possible from all tracers in the void. 
Unsurprisingly therefore, sufficiently small spheres around 
the circumcentre contain no particles at all and n{r) ~ 0 for 
voids of all sizes, even though rimin values are never so small 
and vary with void size. The Voronoi tessellation avoids this 
issue because of its self-adaptive resolution. It can be seen 
from Fig. 5 that the Voronoi reconstructed rimin is a much 
better predictor of the true dark matter density than num¬ 
ber counts — which is why it is preferable to reconstruct 
the density field from the tessellation in the first place. For 
this reason, such stacked number density profiles should not 
be relied upon for quantitative analysis without calibration. 

For this purpose we instead make use of the full dark 
matter density field at high resolution. Profiles of the av¬ 
erage enclosed dark matter density 1 + A(r) are shown in 
Fig. 7, for the same void stacks as before. These confirm 
some properties of watershed voids which are of significance 
for the attempts to model them theoretically. First, as al¬ 
ready seen in Figs. 2 and 3, larger voids contain deeper den- 
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Figure 6. Stacked tracer number density profiles for voids of different sizes. Stacks are chosen to include equal numbers of voids in each. 
Left: Profiles for void stacks centred on void barycentres. The solid lines show the best-fit forms of the fitting formula of Hamaus et al. 
(2014); Sutter et al. (2014) (eqs. 13, A2 and A3), which generally provides a poor fit to the data. Due to discreteness artefacts at small 
r, data points shown with open symbols are excluded from the fitting procedure. Right: Profiles for the same voids but with the stacks 
centred around the circumcentres. By construction the circumcentre more accurately locates the region of minimum tracer density within 
the void. 


sity minima. Secondly, the enclosed density contrast within 
these voids is A(r) > —0.8, for all void sizes and at all dis¬ 
tances r. The condition for shell crossing to occur is thus 
not satisfied at any point within the average void. Finally, 
the central matter densities are much lower for circumcen¬ 
tre stacks than those centred on the barycentre, as expected 
from Fig. 3. 

Our results in Fig. 6 may also be compared to those of 
Hamaus et al. (2014), who postulate a ‘universal’ profile for 
voids based on the functional form 


n(r) 


— 1 -f (5c 


( 1 - jr/r.r \ 

\l + ir/R,f) ■ 


(13) 


This profile form has two free parameters, (5c and r^, with 
a{rs) and /?(rs) fixed by eqs. A2 and A3. Note that although 
Hamaus et al. (2014) refer to this as the density profile p{r), 
their fits to data are in fact based on measurements of the 
tracer number density n(r) — as emphasized above, for sub¬ 
sampled dark matter tracers these two quantities are not the 
same. We determine the best-fit values of these two parame¬ 
ters by fitting to the n(r) data for each void stack. To avoid 
the discreteness artefacts described above, in each stack we 
exclude data points with (r/iZ„) Rv < (3/47r)^^^ rjv from the 
fitting procedure, where Rv is the mean void radius for the 
stack. The resulting fits are shown by the solid lines in the 
left-hand panel of Fig. 6. It can be seen that this ‘universal’ 
profile function generally provides a poor fit to the data, 
both within the void interior and in the overdense walls. 

Similar results are obtained when fitting to the true 
dark matter density profiles p(r). The behaviour of the best- 
fit parameters in equation 13 as functions of the void size 
Rv also significantly differs from that claimed by Hamaus 
et al. (2014) and Sutter et al. (2014). In particular, these 
authors suggest that the central density contrast 5c is an 
approximately linear, increasing function of Rv for voids of 
all sizes and in all tracer populations. This behaviour is cen¬ 
tral to their claim that equation 13 can provide a ‘universal’ 


description of void density profiles in all tracer populations 
through a simple rescaling of void sizes. On the contrary, 
as shown in Fig. 3, 5c{Rv) is a strictly decreasing function 
for Main sample voids, and is a non-linear, U-shaped func¬ 
tion for Dense sample voids. Further discussion of the fitting 
profile is provided in Appendix A, where we show that the 
behaviour of other parameter fits also disagrees with the 
claimed universality. 

We conclude that the fits provided by Hamaus et al. 
(2014); Sutter et al. (2014) fail both qualitatively and 
quantitatively to describe the density profiles we observe. 
Equally, we find no evidence for the self-similarity of tracer 
density profiles seen by Nadathur et al. (2015), but in this 
case differences in methodology, the use of dark matter par¬ 
ticles instead of galaxies as tracers and the different selection 
criteria applied to voids may preclude direct comparison (see 
Nadathur & Hotchkiss 2015 for a fuller discussion). 

So far, following earlier works (Hamaus et al. 2014; Na¬ 
dathur et al. 2015) we have only considered the variation in 
the mean profile with the size of the voids included in the 
stack, but it is clear that this cannot be the only important 
variable. In fact, as shown in Fig. 8, voids of similar sizes but 
different minimum densities rimin have very different density 
profiles. Voids with different rimin clearly do not enclose the 
same density contrasts, and deeper density minima do not 
correspond to steeper density profiles. The enclosed density 
contrast A clearly varies widely over the void population, 
providing further evidence that the population of watershed 
voids does not satisfy the foundational assumption of the 
excursion set model. 

It is also clear that a more complete description of the 
density profiles around voids is obtained by accounting for 
the extent of variation in both dimensions of the (rimin, Rv) 
plane. Fitting formulae such as those provided by Hamaus 
et al. (2014) or Nadathur et al. (2015), which account only 
for variation with void radius, will in principle be unable to 
describe the full variety of watershed voids. 
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Figure 7. Stacked profiles of the total enclosed dark matter density, 1 + A(r), within radius r of the void centre. The stacks are the 
same as in Fig. 6. The left panel shows profiles for stacks centred around the void barycentres, and the right panel for stacks centred 
around the circumcentres. 
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Figure 8. Stacked profiles of the total enclosed dark matter density, 1 +A(r), for voids within the same size range 15 < Rv < 20 h~^Mpc, 
but with different minimum tracer densities. Profiles in the left panel are stacked about the barycentres, and in the right panel about 
the circumcentres. 


However, it is worth stressing that the distribution 
of highly biased galaxies trace dark matter underdensities 
rather differently than the randomly down-sampled dark 
matter particles we have used in this work (Nadathur & 
Hotchkiss 2015), and it is the dark matter prohles of galaxy 
voids which are of greater practical interest in cosmology. 


5 CONCLUSIONS 

Our aim in this paper was to provide an empirical investiga¬ 
tion into the properties of watershed voids in order to better 
understand the operation of void finding algorithms snch as 
VIDE and ZOBOV and the relation to theory. Several pre- 
vions studies have focused on the distribution of void sizes 
alone, and have attempted to ht this using modifications 
of the spherical evolution model. Such an approach however 
misses the important relationship between void size and den¬ 


sity: larger voids correspond to deeper density minima. This 
is a fundamental feature of ZOBOV that holds irrespective 
of whether the tracers used for void identification are simu¬ 
lation dark matter particles, haloes or galaxies. It is also a 
more general property that should apply to any watershed 
void finder. 

The conclusion that follows from this relationship — 
and which we also demonstrate directly through stacked 
density profiles around void centres — is that watershed 
voids cannot correspond to a population of objects which 
all enclose the same density contrast, which is the princi¬ 
pal starting assumption of theoretical descriptions deriving 
from the model of Sheth & van de Weygaert (2004). It has 
long been known that the void number function prediction 
of this model fails to match that of watershed voids by many 
orders of magnitude. It has sometimes been argued without 
proof (e.g. Sutter et al. 2014; Chan et al. 2014) that the void 
formation threshold 5^ might differ from the shell crossing 
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value in the spherical model due to the generally aspherical 
nature of watershed voids. This assumption has led several 
authors to treat 5v as a free parameter but without altering 
the basic model. However, given the range in enclosed den¬ 
sity contrasts A over the watershed void population, a single 
value of (5v for all voids does not seem tenable. Even more 
suggestive is the fact that for no subset of these voids does 
the average enclosed density contrast satisfy the criterion 
for shell crossing, A ~ —0.8, at any radial distance from 
the centre, let alone at the void radius R^. Nor do smaller 
voids correspond to deeper density minima as expected in 
the model. 

The simplest interpretation of this evidence is that wa¬ 
tershed voids simply do not correspond to objects that have 
undergone shell crossing. With hindsight this should not 
seem surprising — ZOBOV uses only information on the lo¬ 
cal topology of the density field, and makes no reference to 
shell crossing. Furthermore, neither VIDE nor ZOBOV ap¬ 
ply any meaningful conditions even on the minimum tracer 
density rimin within voids, instead reporting all local density 
minima. Attempts to explain how the shell crossing density 
criterion may be altered in such voids therefore seem to be 
misguided. A simpler starting proposition would be to give 
up the enforced assumption of shell crossing and to describe 
watershed voids simply as what they are: regions of density 
minima. 

We should stress that breaking this link to theoretical 
models of shell-crossed voids does not necessarily make the 
results obtained from watershed void finders less useful for 
practical cosmological studies. For instance, these voids can 
still be used to identify large-scale underdense environments. 
Some of them (though not all) will also correspond to max¬ 
ima of the gravitational potential, and so they can still be 
used for studies of lensing (Melchior et al. 2014) or the ISW 
effect (Cai et al. 2014; Hotchkiss et al. 2015; Planck Collab¬ 
oration et al. 2015). Equally, we do not intend to claim that 
the Sheth & van de Weygaert (2004) model does not cor¬ 
rectly describe shell-crossed underdensities on much smaller 
scales (although see Falck & Neyrinck 2015; Achitouv et al. 
2015). Our statement is simply that this and related models 
do not match simulation or observational data because the 
word ‘void’ has a different meaning in the two contexts. 

Another interesting feature of our results is the rela¬ 
tionship between the underdensity in voids measured using 
subsampled tracers and using the full resolution dark matter 
density. We show that values of n and p do not completely 
agree, and apparent tracer underdensities in deep voids are 
deeper than the true dark matter minima. The relationship 
between n and p depends on the mean sampling density 
of tracers; it will certainly also change if biased tracers are 
used. This does not affect the basic operation of ZOBOV, 
which only uses relative tracer densities to identify minima, 
but it argues against the use of absolute values of the central 
tracer density in applying selection cuts, as has sometimes 
been suggested (e.g. Sutter et al. 2012; Jennings et al. 2013; 
Sutter et al. 2014; Nadathur & Hotchkiss 2014). In other 
words, selecting a region which apparently satisfies the shell¬ 
crossing criteria in terms of the tracer number density does 
not ensure that it does so in the true matter density. This 
was already pointed out by Furlanetto & Piran (2006) for 
the case when the tracers are galaxies; our results show that 


it applies even if the tracers are a subset of dark matter 
particles in the simulation. 
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APPENDIX A: VOID DENSITY PROFILES 


In this section we provide a quantitative discussion of the 
differences between our results on the void density profiles 
and those of Hamaus et al. (2014); Sutter et al. (2014). For 
this purpose we consider only the average profiles for stacks 
centred on the void barycentres and we also select the 
stacks on the basis of void size Rv alone. Note that this last 
step ignores the very strong systematic dependence of pro¬ 
files on rimin seen in Fig. 8 and by Nadathur & Hotchkiss 
(2015), so is intended solely in order to enable a direct com¬ 
parison with the earlier results. 

We consider both stacked tracer number density and 
matter density profiles for our Main sample voids, which we 
will characterize in terms of the general fitting function of 
Hamaus et al. (2014), 


d{r) 

d 


1 + 


( 1 - {r/rsT \ 


(Al) 


where d = n, p respectively, and the values of the parameters 
Sc, Ts, a and /J are fit to the data in each case. Note that 
Hamaus et al. (2014); Sutter et al. (2014) only provide fits to 
n(r), but refer to this as the true density profile. The form 
of equation Al considered in these works has only two free 
parameters. Sc and r^, since a and P satisfy 


a{rs) ~ -2.0{ra/Rv) + 4.0 


(A2) 


and 

, 17.5(r,/R„) - 6.5 for r./R. < 0.91, , . 

\ -9.8(r,/Rc) -h 18.4 for r,/R„ > 0.91. ^ 

In Section 4.3 we noted that this restricted two-parameter 
form does not in general provide a good fit to the n(r) pro¬ 
files of Main sample voids. The same is also generally true 
for fits to p{r). We therefore treat a and /3 as free parame¬ 
ters and compare their best-fit values thus obtained to those 
predicted by eqs. A2 and A3. 

Fig. Al shows the best-fit values of a and P as func¬ 
tions of Ts/R^ for stacks of different mean void radius R„. 
For comparison, the black dashed lines in each case show the 
dependence expected from eqs. A2 and A3 respectively. The 
stacks are created by binning the void population on the ba¬ 
sis of the void size; the bins are chosen in order to have equal 
numbers of voids and thus equal statistical power to con¬ 
strain the profile. Fits are shown for both the tracer number 
density (green circles) and matter density (red squares) pro¬ 
files, using the appropriate version of equation Al in each 
case. For tracer number density profiles n(r), data points 
with {r/Rc)Rv < {3/4nf^^rN are not considered in ob¬ 
taining the fits, due to the artefacts discussed in Section 
4.3. Fits to the p{r) and n(r) profiles differ significantly from 
each other, as we have repeatedly emphasized, despite the 
tracers being a subset of dark matter particles. Both a and 
P are also clearly discrepant with eqs. A2 and A3, mean¬ 
ing that the two-parameter fitting formula of Hamaus et al. 
(2014); Sutter et al. (2014) is not a good description of the 
observed void density profiles. 

Fig. A2 shows the dependence of the fit values of cen¬ 
tral density contrast he and scale radius Vs on the stack mean 
void radius R„ for our voids. Unsurprisingly n(r) and p(r) 
fits again give very different results. Nevertheless, both sets 
of data points show a clear decrease in Sc with increasing 
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Figure Al. The dependence of the best-fit parameters a and 
/3 in equation Al on the value of for fits to the measured 
tracer number density profiles n(r) (green circles) and true mat¬ 
ter density profiles p(r) (red squares) for stacks of Main sample 
voids centred about the void barycentres. Stacks are created by 
binning the voids according the and bins are chosen to con¬ 
tain equal numbers of voids. The black dashed lines represent the 
expectations if eqs. A2 and A3 were correct. 



Figure A2. The dependence of the best-fit parameters 5^ and 
in equation Al on the stack mean void radius R^, for the same 
stacked data as in Fig. Al. 

Rv The same trend is also seen for voids in galaxy popula¬ 
tions (Nadathur & Hotchkiss 2015). This is in clear disagree¬ 
ment with the results of Hamaus et al. (2014) and Sutter 
et al. (2014), who suggest 5c{Rv) is an increasing function 
for voids in all tracer populations. 

Sutter et al. (2014) also make the stronger claim of uni¬ 
versality: that the fitted values of Sc{Rv) for voids in any 



Figure Bl. The contamination fraction of voids in the Main 
sample. Voids are regarded as genuine if the matter density at 
their circumcentre locations is less than average, and spurious 
if not. The Poisson likelihood contours, minimum size resolution 
and mean inter-particle separation are represented as in Fig. 2. 

tracer population can be translated to the corresponding 
values for voids in any other tracer population, irrespective 
of the sampling density and tracer type, by a simple rescaling 
of the void radius. Such a rescaling is however only possible 
if 5c{Rv) is approximately linear, and shows the same be¬ 
haviour (either increasing or decreasing with for voids 
in all samples. Fig. 3 shows that in the Dense sample, the 
density at the void barycentre is far from a linear function 
of Ry, with a clear minimum at ~ 15 /i“^Mpc. Com¬ 
paring this to Fig. A2 shows that it is not possible to relate 
void profiles from these two samples by a rescaling of the 
void radius. We conclude that the void density profile is not 
universal, but depends at least on the mean density of the 
tracer population. In Nadathur & Hotchkiss (2015) we ar¬ 
gue that it also depends on the tracer bias, and the choice 
of arbitrary input parameters in the watershed algorithm. 


APPENDIX B: EFFECTS OF TRACER SHOT 
NOISE ON VOIDS 

Voids are by definition regions containing a smaller than 
average number of tracer particles, so their identification 
is necessarily affected by shot noise problems. To a very 
large extent, the use of the VTFE reconstruction in the 
void-finding algorithm mitigates this problem compared to 
naive methods of estimating tracer number densities. Nev¬ 
ertheless, there is a real possibility that some fraction of the 
voids returned by the algorithm are spurious detections that 
are not related to true matter underdensities, especially as 
ZOBOV detects large numbers of ‘voids’ even in completely 
random Poisson point distributions (Neyrinck 2008). This is¬ 
sue is particularly important when using sub-sampled dark 
matter tracers. 

To quantify the level of possible Poisson contamination 
we examine the true matter density in the simulation out¬ 
put at the location of the void circumcentre which is 
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the best estimate of the location of minimum density. Voids 
are then classified as ‘spurious’ if they are overdense at this 
location, pcc/p > 1, or ‘genuine’ if pcc/p < 1- This classi- 
hcation deliberately avoids more stringent canonical limits 
on pec based on comparison with the excursion set model, 
since this has been shown not to apply to watershed voids 
in any case. It also uses the circumcentre rather than the 
barycentre, since the barycentre has been shown to be a 
worse locator of the true matter underdensity. 

Based on this classification, we calculate the contamina¬ 
tion fraction. This is shown in Fig. B1 for the Main sample 
as a function of void observables Umin and The con¬ 
tamination understandably increases slightly at larger rimin. 
There is also possibly a small increase in the contamina¬ 
tion fraction within the overlap region with the Poisson con¬ 
tours. However, contamination is generally low and for al¬ 
most all (rimin, Rv) values the majority of voids are genuine. 
It is also clear that P (Poisson| (rimm, is significantly 
smaller than P ((rimin, Pn)| Poisson) within much of the over¬ 
lap region. Comparison of the contamination fractions for 
the Dense and Main samples shows that, as expected, the 
contamination fraction increases with subsampling of the 
tracer population. Nadathur & Hotchkiss (2015) find that 
at the same sampling density contamination decreases with 
increased tracer bias. 

Shot noise due to subsampling of tracers also reduces 
the number of detected voids at almost all R^, as shown in 
Fig. 1. For the voids that are detected, it can also cause an 
offset in the recovered locations of the void centre. To in¬ 
vestigate the relative stability of the alternative centre def¬ 
initions and XT, we matched each void centre in the 
Main sample to its nearest neighbour in the Dense sample 
and examined the distribution of matched nearest neigh¬ 
bour distances in both cases. Subsampling can lead to a 
signiheant offset: for the circumcentre the median of this 
distribution was 0.44P„ and the mean 0.51i?i,. The corre¬ 
sponding values for the matched barycentre locations were 
7 and 3 per cent larger respectively, indicating that the cir¬ 
cumcentre definition is somewhat more robust to shot noise 
effects. We also repeat the matching for voids in different re¬ 
alizations of the same random subsampling: in this case the 
recovered voids can correspond to rather different subsets 
of the original voids, so close matches are less common for 
both centre definitions. However the median matched dis¬ 
tance for circumcentres remained smaller than or equal to 
that for barycentres in all our tests. 
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